/*
 * EBFV1_calculateVelocityField.cpp
 *
 *  Created on: 14/09/2012
 *      Author: rogsoares
 */

#include "EBFV1_hyperbolic.h"

namespace PRS{

// computes total velocity on edges for domain dom
	double EBFV1_hyperbolic::calculateVelocityField(int dom, int dim, TMesh* pTMesh_data){
		int nedges,n;
		// get absolute permeability tensor: dim x dim
		int flag_dom = pTMesh_data->getdomFlags(dom);
		const double *K = pSimPar->getPermeability(flag_dom);

//		ofstream fid;
//		fid.open("PRS.txt");

		pTMesh_data->edge_list[dom].getSize(nedges,n);
		for (int i=0; i<nedges; i++){
			EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);				// get next domain's edge
			VertexData* vdata_0 = edata->vertex_0;
			VertexData* vdata_1 = edata->vertex_1;
			double* grad_0 = vdata_0->grad_p;
			double* grad_1 = vdata_1->grad_p;
			double grad_mean[3];
			for (int j=0; j<dim; j++){
				grad_mean[j] = .5*(grad_0[j] + grad_1[j]);
			}
			double inner = dot(edata->versor,grad_mean,dim);
			double VIJN[3] = {.0,.0,.0};
			for (int j=0; j<dim; j++){
				VIJN[j] = grad_mean[j] - inner*edata->versor[j];
			}

			// velocity parallel component - finite difference approximation
			double v_FD = (vdata_1->p - vdata_0->p)/edata->length;
			//double VIJFD[3] = {-edata->versor[0]*v_FD, -edata->versor[1]*v_FD, -edata->versor[2]*v_FD};
			double VIJFD[3] = {edata->versor[0]*v_FD, edata->versor[1]*v_FD, edata->versor[2]*v_FD};

			// total velocity: vIJ
			double v_total[3] = {VIJFD[0]+VIJN[0], VIJFD[1]+VIJN[1], VIJFD[2]+VIJN[2]};

			// make v_old = v_new (for t=0, v_new = 0)
			for (int j=0; j<dim; j++){
				edata->vel_old[j] = edata->vel_new[j];
			}

			// get mobility for nodes I and J
			double MobIJ = 0.5*(vdata_0->lambda_total + vdata_1->lambda_total);			// average total mobility between vertex 0 and vertex 1

			int pos = 0;
			double vel[3] = {.0,.0,.0};
			if ( pSimPar->is_K_Isotropic() ){
				for (int j=0; j<dim; j++){
					vel[j] = K[pos]*v_total[j];
					pos += dim+1;
				}
			}
			else{
				for (int k=0; k<dim; k++){
					for (int j=0; j<dim; j++){
						vel[k] += K[dim*k+j]*v_total[j];
					}
				}
			}
			for (int j=0; j<dim; j++){
				vel[j] *= -MobIJ;
				if (fabs(vel[j])>1e8){
					throw Exception(__LINE__,__FILE__,"opa!");
				}
				edata->vel_new[j] = vel[j];
				edata->vertex_0->vel[j] += vel[j];
				edata->vertex_1->vel[j] += vel[j];
			}
			//char text[256];
			//sprintf(text,"vel[%d - %d]: %.8E\t%.8E\tgrad_mean: %.8E\t%.8E\n",vdata_0->ID,vdata_1->ID,edata->vel_new[0],edata->vel_new[1],grad_mean[0],grad_mean[1]);
			//sprintf(text,"vel[%d - %d]: VIJFD: %.8E %.8E\n",vdata_0->ID,vdata_1->ID,VIJFD[0],VIJFD[1]);
			//fid << text;
		}
//		fid.close();
//		STOP();
	}
}
